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ABSTRACT 

We describe a Monte Carlo radiative transport code intended for calculat- 
ing spectra of hot, optically thin plasmas in full general relativity. The version 
we describe here is designed to model hot accretion flows in the Kerr metric 
and therefore incorporates synchrotron emission and absorption, and Compton 
scattering. The code can be readily generalized, however, to account for other 
radiative processes and an arbitrary spacetime. We describe a suite of test prob- 
lems, and demonstrate the expected N^ 1 ^ 2 convergence rate, where N is the 
number of Monte Carlo samples. Finally we illustrate the capabilities of the 
code with a model calculation, a spectrum of the slowly accreting black hole Sgr 
A* based on data provided by a numerical general relativistic MHD model of the 
accreting plasma. 

Subject headings: numerical methods, radiative transfer, magnetohydrodynamics 

1. Introduction 

There is wide interest in calculating the emergent radiation from relativistic astrophys- 
ical sources, including accreting black holes, accreting neutron stars, and relativistic blast 
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waves. A variety of methods for solving the r adiative transfer problem in these sources have 



been developed over the last few c 


ecades (e.g. 


Pozdvnakov et al. 1983 


Gorecki & Wilczewsk] 




1984; 


Hauschildt & Wehrsd 


1991; 


Carrigan & KatJl992: CoDDi et al. 


1993; Stern et a 


1995 


Poutanen & Svensson 


L996: 


Zane et al.l 1996; 


Dove et al.lll997l: Bottcher & LianekoOl: 


Schnittman et al. 


20061; 


Noble et al. 2007 


; Wu et al. 


2008), some based on Monte Carlo schemes. Few of these 



schemes take full account of relativistic effects in the source, however, and this is crucial in 
estimating the spectra of hot plasma deep in a gravitational potential well, or highly rela- 
tivistic blast waves. Monte Carlo t ransport of radiation in accretion flows around compact 



objects has been considered by (e.g. ISchnittman &; Krolik 



2005 



1996 



Stern et al 



Bottcher et al.ll2003l ; IBottcher fc Liang 



2001 



2009 



Schnittman 



Laurent fc Titarchuk 



2006 



1999 



lYao et al. 



Agol &; Blaes 



19951). Amo ng oth ers. ICullenl d200lh: iMolnar fc Birkinshawi Jl999h: iHua 



19971 ) ; I Gorecki h Wilczewskil (119841 ) ; iPozdynakov et al.l (119831 ) give more general discussions 
of Monte Carlo radiative transfer techniques. 

We were motivated by efforts to model the radio source and black hole candidate Sgr A*. 
Our interest in this source drove us to develop a numerical scheme that could accurately cal- 
culate spectra of a relativistic source in which the plasma properties (velocity, density, mag- 
netic field strength, and temperature) were specified by a separate model — that is, sources 
in which radiation plays a negligible role in the dynamics and energetics. The result, a 
Monte Carlo scheme called grmonty, is described in this paper. The spirit of our calculation 
is to obtain an accurate spectrum with as few approximations as possible. To this end we 
treat Compton scattering with no expansions in v/c, and allow for general angle-dependent 
emission and absorption (we specialize to thermal synchrotron in this work). 

In designing grmonty our philosophy has been to maximize the physical transparency 
and minimize the length of the code, occasionally at the cost of reduced performance. Some- 
times simplicity and efficiency are in harmony. We chose to directly integrate the geodesic 
equation rather than using a scheme that relies on integrability of geodesies in the Kerr 
metric. We will show that for radiative transfer problems where many points are required 
along each geodesic, direct integration is not only simpler and easier to modify, but also 
faster. 

Our paper is organized as follows. In §2] we describe how we sample emission, and in 
§|3] we describe how we track photons along geodesies. Evolution of superphoton weights 
under absorption is described in §HJ and sampling of scattered photons is discussed in §[5j 
Photons at large distance from the source must be sampled and assembled into spectra; 
this is described in §EJ The code has been extensively verified; we describe tests in §7J §S] 
describes a sample calculation, and $9] summarizes our results. 



Throughout this paper we assume that there is an underlying model that can be queried 
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to supply the rest-mass density po, the internal energy u, the four-velocity u^, and magnetic 
field four-vector for the radiating plasma. Usually we expect the model to be supplied by 
a numerical simulation in a coordinate basis x M . 



2. Manufacturing Superphotons 

Emission in grmonty is treated by sampling the emitted photon field. The samples, 
here called "superphotons" (also "photon packets"), have weight w, coordinates and 
wave vector The weight w ^> 1 is a pure number that represents the ratio of photons to 
superphotons: dN = wdN s (N s = number of superphotons, N = number of photons). In our 
models the weight is a function of the emitting plasma frame frequency v and nothing else. 
The coordinates are typically in model units (e.g. for a black hole accretion flow calculation, 
length unit L = GM/c 2 ), and the components of k^ are given in units of m e c 2 . 

How should superphotons be distributed over x^ and k^l The initial superphoton 
momentum can be described in an orthonormal tetrad basis >, that is attached to the 

(a) 

plasma, so that = uF" (fi is the coordinate index, and (a) is the index associated with 
the tetrad basis, raised and lowered using the Minkowski metric). In the tetrad basis k^ 
is specified by frequency v and spatial direction unit vector n that is contained within the 
solid angle dfl. The probability distribution for superphotons is then 

1 dN s 1 dN 1 j v 



-g d 3 xdt dvdVL w^f^g d 3 xdt dvdVL w hv 



(1) 



where j v is the emissivity (always defined in the plasma frame), since \f—g d 3 xdt is invariant 
(meaning coordinate invariant). In a time interval At we expect to create 

N sMt = At ! d 3 x dv dQ-^- (2) 

J w hv 

superphotons over the entire model volume. The total computational effort is proportional 
to N s>tot , so we control the computational effort by scaling the weights. 

How should we distribute superphotons over the volume? grmonty subdivides the model 
volume into volume elements ("zones") of size A 3 x (e.g. in Boyer-Lindquist t,r,8,<f) coordi- 
nates for the Kerr metric, ArA9A<f)). For zone i we expect 

N s i = At A 3 rrv^ / dvdQ-^- (3) 
J w hv 

to be created in time At. We create N Sji superphotons at the center of zone i. Fractional 
N s ^i are dealt with by rejection sampling. 
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The momentum-space (wave vector) piece of the probability distribution ([I]) can be 
sampled by techniques outlined below to give v and n. With x^, v and n in hand, we can 
construct k^ and, finally, transform to the coordinate basis efsk^ = k^. 

The remaining ingredients in the sampling procedure are the emissivity, the orthonormal 
tetrads, and the sampling procedures for v and n. 



2.1. Emissivity 



grmonty depends on the emissivity only through functions that specify j v and J dudQj u / u, 
so it is straightforward to include any emission/absorption process. 

In our target problem the only source of su perphotons i s therm al synchrotron emission 
at dimensionless temperature e = kT e /(m e c 2 ). iLeung et al.l (120091 ) show that, for e > 0.5, 



U- 



eB 



9 \27im e c 



(4a) 
(4b) 
(4c) 



where K 2 is the modified Bessel function of the second kind, n e is the number density of 
electrons, B is the magnetic field strength, and 9 is the angle between the wave vector and 
magnetic field . For large 9 P , Kpj Oj 1 ) ~ 26^, but for 6 e < 1 better agreement with the 
emissivities of ILeung et al.l (120091 ) is obtained if K 2 is evaluated directly. Since the emissivity 
must be evaluated many times, it is most efficient to precompute ^(©eT 1 ) at the beginning 
of the calculation and store the results in a table. 



2.2. Orthonormal tetrads 

The wave vector sampling is done in an orthonormal tetrad attached to the fluid. 
We construct the orthonormal tetrad using numerical Gram-Schmidt orthogonalization. 
Here \x is the coordinate index, and (a) is the index associated with the tetrad basis, which 
is raised and lowered using the Minkowski metric. 

We set e/L = w M {u^ = plasma four- velocity), and then use 6 M , the magnetic field four- 
vector, as the first trial vector (this is numerically convenient since we will want to orient 
wave vectors with respect to the magnetic field; if = then we use a default, radius- 
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aligned, trial vector). Thus efa = NORM [b^ — efoiefah)) (NORM: normalize). The 
process is repeated with additional trial vectors to create a full tetrad basis. 

The tetrad-to-coordinate basis transformation is 

k? = ef a) *M (5) 
and the coordinate to tetrad transformation is 

jfeW = eWjfe". ( 6 ) 

With these transformations in hand, we can construct the superphoton wave vector in the 
orthonormal frame and then transform it to the coordinate basis. 



2.3. Wave vector sampling procedure 

2.3.1. Photon energy 
Within zone i superphotons are distributed over frequency according to the distribution 

f/As '' AtA 3 x^-^- [dQj u . (7) 



d In v hw 



max ■ 



We sample the distribution only between a minimum and maximum frequency v m i n and v. 
These must be chosen so that no significant emission is omitted from the final spectrum. 

We distribute superphotons over frequency by rejection sampling. For simplicity, we 
use a constant envelope function equal to the maximum of Equation (j7j) (for zone i). Thus 
we draw tentative values uniformly in Inu from v m i n to v max . The efficiency of the sampling 
procedure is the ratio of the areas under the distribution and envelope, so if the distribution 
given in Equation §J§ is sharply peaked this technique can be inefficient. 

In practice, we choose a tentative frequency v = exp(rx In v max j u min + \n. v min ), where r 1 
is drawn from a uniform distribution on [0,1) (we use the Mersenne twister random number 
generator from the GNU Scientific Library, hereafter GSL). A second number T2 is drawn 
from [0,1) and the process is repeated until 



dN 

r 2 < 



dlnv 



MAX 



dN s 
din 



v 



The efficiency of this process is ~ 15%, but the cost is small compared to the total cost of 
grmonty. 
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2.3.2. Photon direction 

The superphoton direction n is described by polar coordinates 9 and <fi in the tetrad 
frame, where is the angle between the spatial part of the wave vector and the magnetic 
field. The colatitude 9 is obtained by rejection sampling: a tentative value for 9 is obtained 
by drawing \i = cos9 from a uniform distribution on [-1,1), a second number r is drawn from 
a uniform distribution on [0,1), and 9 is accepted if 

Ju(n/2) 1 ; 

(this procedure is specific to the synchrotron emissivity). The efficiency of this scheme is 
problem dependent; for our target application the efficiency is ~ 65%. Finally, <f> is drawn 
from a uniform distribution on [0,27r). 



2.3.3. Transformation to coordinate frame 

Once 9, 0, and e = hu/m e c 2 are selected, the wave vector is completely specified in the 
orthonormal tetrad frame: 

A; (0) = e (10a) 
k( 1 ) = eC os9 (10b) 
k^ = esin^cos^ (10c) 
k^ = esin6>sin0, (lOd) 

and the wave vector in the coordinate frame is = ef ,k^ a \ 

(a) 



3. Geodesic Integration 

General relativistic radiative transfer differs from conventional radiative transfer in 
Minkowski space in that photon trajectories are no longer trivial; photons move along 
geodesies. Tracking geodesies is a significant computational expense in grmonty. 

The governing equations for a photon trajectory are 

which defines A, the affine parameter, the geodesic equation 

dk a 



= ( 12 ) 



d\ 
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and the definition of the connection coefficients 



(13) 



in a coordinate basis. 



We assume nothing about the metric, so it is easy to change coordinate systems and 
even to extend the code to dynamical spacetimes. Nevertheless, our main application — to 
black hole accretion flows — is in the Kerr metric, where geodesies are integrable. The four 
constants of the motion are the energy- at- infinity E (in Boyer-Lindquist coordinates t, r, 9, 4>, 
E = —kt), the an gular momentum I = k^, Carter's constant Q = kg + k\ cot 2 9 — a 2 k 2 cos 2 9 



sec 



Carterlll968l ). and the condition that k a be null: k^k^ = 0, equivalent to the dispersion 
relation for photons in vacuo: uj 2 = c 2 k 2 . These four constants of the motion can be used to 
quasi-ana lytically obtain x^ and k^ in terms of an initial (or final) position and w ave vector 
(see, e.g., iRauch & Blandfordl[l994l : beckwith fc Done!l2005l : IPexter & Agoll2009h . 



The integrability of geodesies in the Kerr metric would appear to provide an opportunity 
for significant computational economies. We show below, however, that direct integration of 
equations (fTTj) and (fT2"|) is not only simpler and more flexible but also faster than at least 
one implementation of an integral-based technique. 

Which ODE integration algorithm is best for the geodesic equation? If only a few coor- 
dinate evaluations are required over the entire geodesic then a high order scheme is optimal. 
For example, we have found that the embedded Runge-Kutta Prince-Dorman method avail- 
able in GSL is fast and accurate; it can easily be made to conserve the integrals of motion to 
machine precision. Many coordinate evaluations are required, however, when integrating the 
equation of radiative transfer, as grmonty does, along superphoton trajectories. A second 
order scheme can then provide the required accuracy at minimal cost. 

Evaluating the connection coefficients is expensive, so we want to choose a scheme that 
minimizes the number of evaluations. The velocity Verlet algorithm, which for the geodesic 
equation is 



x 



n+l 



AA 



dk a 



- r a „ , (x n+ 1 ) k^ + 1 k u 



+l,p n+l,p 



n+l 



h a - b a A 



dk a 



+ 



dk a 



n+l 



(AA), 



(14a) 
(14b) 
(14c) 
(14d) 
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requires only one evaluation of the connection coefficients per step. Accuracy can be im- 
proved by using the result of equation flf 4dj) to recompute the derivative equation (I14cl) 
with k% +lp = k^ + i and then reevaluating equation (114dj) . This process is repeated until 
the change in the wave vector between estimates is below some tolerance. In grmonty we 
continue this iteration until the fractional change is less than 10 -3 (typically only once or 
twice). This does not require any additional evaluations of r° . Very rarely we find that 
this iteration fails to converge, and then grmonty defaults to taking the step with a classical 
4 th -order Runge-Kutta technique. 

How fast and accurate is our geodesic integration scheme? We propose the following 
benchmark. Consider a point on a direct, circular, marginally stable orbit in the equatorial 
plane of a black hole with spin a/M = 1— 2 -4 = 0.9375 that emits radiation isotropically in its 
rest frame. Sample the emitted photons (in a Monte C arlo sense; the analytic circular orbit 



orthonormal tetrads available in iBardeen et al.l (119721 ) are useful for constructing the initial 



wave vectors) and track them until they cross the horizon or reach rc 2 / (GM) = 100 (r is the 
Boyer-Lindquist or Kerr-Schild radial coordinate). Figure [1] shows as dots a representative 
sample of photon geodesies from grmonty in the coordinate frame, illustrating the effects of 
relativistic beaming, lensing, and frame dragging. 

Second order convergence of the velocity Verlet integration scheme is demonstrated in 
Figure [2], which plots the average fractional error in E, I, and Q as a function of a step-size 
parameter e. We typically set e = 0.04 as a compromise between performance and accuracy; 
the average fractional errors at the end of the integrations are ~2x 10 -3 , ~4x 10 -2 , and 
~ 8 x 10 -2 for E, I, and Q, respectively. We have verified that this choice makes geodesic 
tracking errors subdominant in the error budget for the overall spectrum. 

With e = 0.04, grmonty integrates ~ 16, 700 geodesies sec -1 on a single core of an Intel 
Xeon model E5430. If we use 4 th -order Runge-Kutta exclusively so that the error in E, I, 
and Q is ~1000 times smaller, then the speed is ~ 6,200 geodesies sec -1 . If we use the 
Runge-Kutta Prince-Dorman method in GSL with e = 0.04 the fraction error is ~ 10 -10 
and the speed is ~ 1, 700 geodesies sec" 1 . These result s can be compared to the publicly 



available integral-based geokerr code of iDexter fc Agoll (120091 ). whose geodesies are shown 



as the (more accurate) solid lines in Figure [TJ If we use geokerr to sample each geodesic 
the same number of times as grmonty (~ 180), then on the same machine geokerr runs at 
~ 1,000 geodesies sec -1 . It is possible that other implementations of an integral-of-motion 
based geodesic tracker could be faster. 

If only the initial and final states of the photon are required, we find that geokerr 
computes ~ 77,000 geodesies sec -1 . The adaptive Runge-Kutta Cash-Karp integrator in 
GSL computes ~ 34, 500 geodesies sec -1 with fractional error ~ 10 -3 . 
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4. Absorption 

grmonty treats absorption deterministically. We begin with the radiative transfer equa- 
tion written in the covariant form 

1(1 (IA fk)-( uav j(±). (is) 



1 V 

C dX \ v 3 J \v 2 ) v " <y,u/ V i>3 



sec 



Mihalas fc Mihalaslll984l ). Here l v is specific intensity and a u>a is the absorption coef- 
ficient (which is always evaluated in the fluid frame). The absorption coefficient must be 
computed by a separate subroutine; for thermal synchrotron emission we set a v>a = j v / B v . 
C is a constant that depends on the units of fc M (in grmonty, electron rest mass), v (Hertz), 
and the length unit L for the simulation in cgs units. For grmonty 

Lh , , 

C = - 2 . 16 

m e cr 

Each quantity in parentheses in equation ({15]) is invariant; I v j z/ 3 , for example, is proportional 
to the photon phase space density. 

Since l v j u 3 is proportional to the number of photons moving along each ray, I v j z/ 3 oc w, 
and equation fTl5|) implies (ignoring emission) 

JT = ~ W (17) 

where 

dr a = (va£ Uj a) CdX (18) 
is the differential optical depth to absorption and the quantity in parentheses is the "invariant 
opacity." This equation we integrate with second order accuracy 

r a = - ((va u<a ) n + (ya Vt a)n+i) CAX, (19) 

and then set 

w n+1 = w n e- Ta . (20) 

Since the components of are expressed in units of the electron rest-mass energy, v = 
—k^u^m^/h. Storing the invariant opacity at the end of each step saves computations 
since it can be reused as the beginning of the following step. 



5. Scattering 

Our treatment of scattering consists of two parts: the first determines where a super- 
photon should scatter and the second determines the energy and direction of the scattered 
superphoton. 
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5.1. Selection of scattering optical depth 

When a superphoton is created or scattered grmonty selects the scattering optical depth 
t s at which the next scattering event will take place. Scattering follows the cumulative 
probability distribution 

p=l-e- T ° = r s + 0(r 2 ), (21) 

so superphotons will experience on average r s scattering events when t s < 1. In optically 
thin sources this would result in poor signal to noise in portions of the spectrum dominated 
by scattered light. To overcome this, we use the biased probability distribution 

p = l- e- bTs (22) 



where b is a bias parameter. This technique was originally proposed by iKahnl (119501 ) in the 
context of deep penetration of neutrons in radiation shielding and has since been extensively 
explored in the nuclear engineering literature (often refered to as exponential biasing or ex- 
ponential transform). Whereas in deep penetration problems b < 1 in order to allow for 
sampling of radiation at high optical depths, here we set b > 1 in order to better sample 
scattered photons at low optical depths. Superphotons now experience on average br s scat- 
tering events. Two superphotons emerge from a scattering event: the incident superphoton 
of weight w and a new scattered superphoton. For conservative scattering the incident su- 
perphoton has its weight reset to w(l — 1/b) and the new superphoton has weight w/b, so 
that weight (photon number) is conserved. 

What should we choose for the bias parameter 6? The goal is to set b such that scat- 
tering is more likely to occur in regions which contribute most to the spectrum. This is an 
example of the more general technique of importance sampling. Typically we set the bias 
parameter b = MAX(1, aQg/r S)mttX ) (a is a scaling factor we set to l/(6 e ) 2 where (6 e ) is the 
volume averaged dimensionless temperature and T s ^ max is an estimated maximum scattering 
optical depth) to improve sampling on the high energy side of each scattering order, which is 
populated by photons scattered from high temperature plasma. If the bias factor is too large 
a "chain reaction" results in an exponentially growing number of superphotons; l/r s ^ max is 
an estimate of the critical bias factor in a e = 1 plasma. 

We evaluate the scattering optical depth along geodesies in a manner analogous to the 
absorption optical depth; 

T s = - {{va v , 8 )n + (fcOn+O ( 23 ) 



is the scattering depth along a step. 
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5.1.1. Covariant evaluation of extinction coefficient 



In our applications electron scattering dominates. The general, invariant expression for 
the rate of binary interactions dN ab between a population of particles dN a and dNj, is 

dN ab 1 f d 3 p a d 3 p b dN a dN b 



(24) 



y/—g d 3 xdt 1 + 5 a b J V~9Pa V~9Pb d 3 xd 3 p a d 3 xd 3 pb 

where 5 a b prevents double-counting if a = b, d 3 p = dpidp 2 dp 3 , a is the invariant cross 
section, and v ab = c( l + m 2 ml / ( —PapP^) 2 ) 1 / 2 . This is the manifestly covariant generalization 
of equation (12.7) of lLandau fc Lifshitzl (119751 ). 



We want to use this expression to find the cross section for a photon with wavevector 
&o interacting with a population of particles of mass m > 0. We therefore set dN^/d 3 xd 3 p = 
5{k^ — ho) and, dropping the subscripts on k and p, the integral reduces to 

d 3 p dn m 



dN, 



ac 



(25) 



\J—g d 3 xdt J \f—gp l d 3 p k l 

where dn m = dN m /d 3 x. In Minkowski coordinates {y/—g = 1), define (3 m = the particle 
speed in the plasma frame and \x m = the cosine of the angle between the particle momentum 
and photon momentum in the plasma frame. Then 



dn 



m-y 



d 3 p 



dn r 



[1 - Hmfim) ac - 



dt J d 3 p 
It is convenient to rewrite this rate in terms of a "hot cross section" 

Oh = — [ d 3 p ( ^ ! - (1 - fXmPm) cr 

n m J d 6 p 

so that the interaction rate for a single photon is n m OhC and the extinction coefficient is 



(26) 



(27) 



a v = n m a h . 

So far we have assumed nothing about the interaction process. 



(2* 



5.1.2. Electron scattering 

For electron scattering the cross section is the Klein-Nishina total cross section expressed 
in terms of the photon energy in the electron rest frame = e e = ej e (l — /i e /3 e ), (j e = 
(1 — P 2 )^ 1 ^ 2 and we have substituted the subscript e for m): 

--^( 2 + TO + i ^ log(1+2£j )' (29) 
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Here <jt is the Thomson cross section. For e< 1, 

a KN = a T (l - 2e + 0(e 2 )) (30) 

which is numerically stable for small e, unlike equation (1291) . 
Typically we assume a thermal electron distribution, 

dn e n e "fl(3 e ( 7 



d- h e c K 2 (e-^{-0 e )' (31) 

and evaluate (1271) by direct integration to obtain cr^(O e , e). It is efficient to store the resulting 
cross sections in a two-dimensional lookup table at the beginning of the calculation. Our an 
agrees with Wienke (1985). 



5.2. Scattering kernel 

Once it is determined that a superphoton should be scattered at an event the 
superphoton is passed to a scattering kernel which processes the scattering event according 
to the following procedure. Only unpolarized light is considered. 

First, a plasma frame orthonormal tetrad is constructed by the same Gram-Schmidt 
orthogonalization procedure described in and the old superphoton wave vector is trans- 
formed from the coordinate frame to the tetrad frame. 



Second, a scattering electron is selected. We use the procedure described by lCanfield et al 



( 119871 ). which selects the four-momentum of the scat tering electron with an efficiency of 



72% at e e = 1 and nearly 100% for e < 1 or e > 1 flCanfield et al.lll987h . 



Third, we boost from the plasma frame tetrad to an electron frame tetrad q?-. and 
construct the scattered photon wave vector in this frame. The differential scattering cross 
section is sampled for the scattered photon energy e' e and the scattering angle 9. For low 
energy photons (e e < ei; in grmonty q = 10~ 4 ) the scattering is approximately elastic, so 
we set e' e = e e and sample the Thomson differential cross section 

2n ddT 3^ 2 



<tt d cos 9 8 



(1 + cos^) (32) 



for the scattering angle 9 using a rejection scheme. For e e > ei we sample the Klein-Nishina 
differential cross section 

2-k da K N 1 ,e e e' e 2 



a T de' e e\ e e 



+ --l + cos 2 fl), (33) 
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for e' e using a rejection scheme. Here cos 9 = 1 + l/e e — l/e' e . This procedure is inefficient 
for e e ^> 1, but in our target application such large photon energies are rare. Drawing a 
final angle <j) from a uniform distribution on [0, 2n) completes the specification of the photon 
wave vector 

fc(°> = e' e (34a) 
fcM = e ' e cos (34b) 
e^, sin 8 cos (34c) 
A; (3) = e e sin#sin0 (34d) 

in the electron- frame tetrad basis q^ a y, is aligned parallel to the spatial part of the 
incoming photon wave vector. 

Finally, we boost back from the electron-frame tetrad to the plasma frame tetrad (some 
of these steps are combined in our code for computational efficiency), and use the plasma 
frame basis vectors to obtain the coordinate frame scattered photon wave vector k'. 



6. Spectra 

Spectra can be measured using a "detector" with area AA at distance R, frequency 
channels of logarithmic width A In u, and integration times AT. The flux density in frequency 
bin % is then just 

F -%Aln^ATAA^ (/ ^ (35) 

3 

where the sum is taken over all superphotons j that land in the channel during the inte- 
gration. In principle a software detector can behave just like a physical detector, producing 
time-dependent spectra from time-dependent flow models. In practice time-dependent mod- 
els are not (yet) treated self-consistently. 

To see why, consider a time-dependent model based on a general relativistic magneto- 
hydrodynamics (GRMHD) model of a black hole accretion flow. Self-consistent treatment of 
the radiation field would require generating and tracking superphotons through the simula- 
tion data as it evolves, i.e. coupling the grmonty to the GRMHD code (the simulation data 
could be stored and post-processed, but this would require storing almost every timestep and 
would be impractical and inefficient). The mean number of superphotons tracked simultane- 
ously would depend on the desired signal-to-noise in the final spectrum, as well as the time 
and energy resolution. Our experience suggests that for nominal energy resolution, signal- 
to-noise, and time resolution of order the horizon light crossing time, ~ 10 8 superphotons 
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would need to be maintained within the simulation for the duration of the evolution. The 
total number required is then N tot 2> 10 8 and is currently inaccessible without significant 
computational resources. We plan to attempt this calculation later. 

For now we construct spectra of time-dependent data using a stationary-data (or "fast- 
light") approximation: each time-slice of data is treated as if it were stationary (time- 
independent). The slice emits superphotons for a time At. The photons then propagate 
through the slice data (as t varies along a geodesic the fluid variables are held fixed) and are 
detected at large distance. The time-steady flux is obtained by substituting At for AT in 
equation fl35l) . 



6.1. Measuring vL v 

In practice we measured vL vi rather than F v ^. Since vL v = 4itR 2 i>F u , and R 2 /AA is 
the solid angle Afl occupied by the detector, 



An 



AttAt Ainu 

j 



— Y, W M- ( 36 ) 



Typically our "detectors" capture all the superphotons in a large angular bin AQ around 
the source. For example, in studies of axisymmetric black hole accretion flows the angular 
bins capture all photons with Boyer-Lindquist r > lOOGM/c 2 , 9 n < 9 < 9 n+ \, independent 
of <j), where 9 n are the bin boundaries. 

We can also estimate average values for any quantity Q associated with emission from 
a source (e.g. the absorption optical depth). We define the weight-averaged value of Q via 

(Q) = f^, (37) 
where the sum is taken within an energy and angular bin. 



6.2. Optimal weights 

The fractional variance in the Monte Carlo estimate for vL v is proportional to w 2 / (N S W 2 ), 
where overline means expectation value and N s is the number of superphotons in the bin. 
Evidently the optimal weighting of superphotons is achieved when (1) the weights of super- 
photons are the same within bins, and (2) the superphotons are evenly distributed across 
frequency bins. 
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If, as we created new superphotons, we knew vL v , then we could set 

w v = =t- L„Alni/ (38) 
N s h 

and set N s = N S;tat /N b , where N b is the number of frequency bins. 

Of course we do not know L u , but for the special case of an optically thin emitting 
plasma we can estimate it: 

U~ J d 3 x^/=g J dQ j u (39) 

This estimate assumes that all photons escape to infinity, and it ignores Doppler shift, 
gravitational redshift, scattering, and angular structure in vL v . Nevertheless it is useful 
because (1) it can be calculated before the Monte Carlo calculation begins; (2) it is far 
better to use the information contained in this rough estimate of the spectrum than to 
proceed using, e.g., uniform weights. 



7. Tests 

We verify the accuracy of grmonty by comparing spectra produced on idealized problems 
against a reference spectrum (vL u ) re f computed analytically (when possible) or computed 
by an independent code. For all tests we use the following error norm, which effectively 
measures the maximum of the fractional error, compared to the reference solution, over 
frequency: 

(e) = — l — l( ^V nty r (Z/jMr£/l ^ In v (40) 

where Ainu = ln(u max / u m i n ) is the range of integration. The range of integration is the 
same as plotted in the spectra for each test. 



7.1. Optically thin synchrotron sphere 

First we consider emission from a homogeneous, optically thin spherical cloud of unit 
volume threaded by a vertical magnetic field in flat space. The cloud parameters are O e = 
100, B — 1 G, and n e = 10 15 cm~ 3 , which gives an optical depth at v = 10 9 Hz of ~ 10~ 2 
perpendicular to the magnetic field. The emissivity and absorptivity are constant along any 
line of sight, so 

L = ^(l-e^ L )^ 3u L + 0(r 2 a ) (41) 
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where L is the path length through the sphere. We numerically integrate this expression 
over detector solid angle to compute the spectral energy distribution, and compare with 
the spectrum grmonty produces in Figure [3J Evidently the result is unbiased. Figure @] 
demonstrates that grmonty converges on the correct solution as oc iV -1 / 2 . 



7.2. Optically thick synchrotron sphere 

The second test is identical except that the electron number density and therefore the 
optical depth are increased by a factor of 10 5 . The intensity along any line of sight is again 

I v = — (l-e- a » L ) (42) 

OL v 

The emitting region becomes optically thick when (j u /B u )L > 1. For the thermal syn- 
chrotron emissivity we use, this occurs below a critical frequency. 

The spectrum is shown in Figure [5] and the convergence is shown in Figure [6j The 
figures make two key points: convergence is slow for small numbers of superphotons; and 
the overall magnitude of the error is larger than in the optically thin case shown in Figure HI 

The slow initial convergence is due to the large optical depth at some frequencies. 
When the optical depth is large no superphotons of appreciable weight are recorded until 
some superphotons have been created in the fraction ~ 1/r of the volume that lies within 
the photosphere. Our problem has r ~ 10 5 at v ~ 10 8 Hz. Since (e) effectively measures 
the maximum of the error over frequency, it is not surprising that grmonty requires ~ 10 6 
superphotons before it begins to converge as N^ 1 ^ 2 . 



7.3. Comptonization of soft photons in a spherical cloud of plasma 



This test is based on a problem posed in §6 of lPozdynakov et al.l (119831 ): the spectrum of 
a spherical, homogeneous, unmagnetized cloud of radius R that contains thermal electrons at 
density n e and temperature T e that scatter light from a central, thermal source of temperature 
T s . Absorption and emission in the cloud are neglected. The dimensionless parameters of 



the problem are O e = kT e / (m e c 2 



RaTTie, and Q s = kT s /(m e c 2 



For this test our reference spectrum is computed with an implementation of Pozd- 
nyakov et al.'s Monte Carlo scheme kindly provided by S. Davis. This code, sphere, has 
been modified in two wa ys: we have replaced the approximate hot cross sections defined in 
Pozdynakov et al.l (119831 ) with our more exact, numerically integrated values, and we use 
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the exact Klein-Nishina cross section equation f[29|) when choosing the electron with which 
a photon should scatter. Without th ese changes to sphere dif ferences between the spectra 
are < 1%, consistent with the error iPozdynakov et all (119831 ) quote for their approxima- 
tions. These small differences are enough, however, to prevent grmonty from converging as 
expected for large numbers of superphotons. 

Figures d [8] and [9] show the radiation spectra (upper panels) produced by both codes and 
a fractional difference (bottom panels) between them for O e = 4, Q s = 10 -8 and for various 
values of the optical depth r = 10~ 4 , 0.1, and 3. Figure [TD] demonstrates that grmonty 
converges to the reference solution as oc iV" 1 / 2 for each optical depth. 



7.4. Synchrotron self-absorbed spectra in black hole spacetimes 

We consider two idealized problems: (1) smooth, spherically symmetric infall onto a 
Schwarzschild black hole; (2) a snapshot of a turbulent accretion flow around a Kerr black 
hole with a/M = 0.9 375 produced by general relativistic MHD simulation with the HARM code 



(iGammie et al.ll20031). In this test, our reference solution is computed by the ray tracing code 
ibothros (jNoble et al. 2007 ). ibothros solves the invariant form of the radiative transfer 



equation along geodesies that terminate in a fictitious camera at large distance. Spectra 
are constructed by imaging the source at successive frequencies and performing an angular 
integral over the images to estimate vL v . In these tests scattering is turned off in grmonty. 

There are algorithmic differences between grmonty and ibothros that lead to differences 
in their spectra. 

First, grmonty measures the flux in energy and angular bins, whereas ibothros mea- 
sures the flux at a particular inclination and energy. This is not an important effect unless 
there is sharp angular or energy structure in the spectrum. 

The next difference is more subtle and is related to the treatment of gridded model data 
used to construct these tests. In grmonty quantities such as the density, temperature, etc., 
are viewed as the average of these variables over a grid zone. In ibothros the grid variables 
are viewed as zone-centered samples and a continuous distribution is created by multi-linear 
interpolation between zone centers. The difference is illustrated in one dimension in Fig- 
ure [TTJ ibothros and grmonty therefore differ in zone-averaged emissivity by 0(Ax/L) 2 , 
where Ax is the zone size and L is the characteristic scale of the emitting structure. 



Differences in the grmonty and ibothros spectra of structures with Ax <C L are there- 
fore small. High frequency synchrotron emission, however, is exponentially dominated by 
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emission from a few zones with highest v s (see equation (14b j) ; these are the zones with highest 
temperature or strongest magnetic field). Then Ax/L ~ 1 for high frequency emission, and 
the grmonty and ibothros spectra differ by of order unity. A similar effect occurs for low 
frequency synchrotron emission where the optical depth is large. The spectrum is sensitive 
to the run of physical variables through the photosphere, which has size comparable to or 
smaller than a grid zone. The subsequent test results will omit these parts of the spectra. 
These inconsistencies between grmonty and ibothros could be eliminated by using identical, 
continuous models for subgrid reconstruction. But this would require significant investment 
in recoding that, in our view, is not worthwhile: to the extent that the spectrum depends 
on the flow structure at and below the grid scale, it is not reliable! 

In spite of these differences, and other differences between the codes related to differing 
accuracy parameters, grmonty should converge on the ibothros result until the effects of 
data interpolation become the dominant sources of error. 

Figures [121 and [T3l show the spectrum and convergence relative to ibothros, respectively, 
for the spherical accretion problem. This problem is an attractive test for at least two 
reasons. First, the emission is isotropic so the effects of angular binning in grmonty are 
eliminated. Second, the flow is smooth, so the differences due to data interpolation will be 
small. Evidently the spherical accretion model converges at the expected, iV -1 / 2 , rate. 

Figures [T41 and [TBI show the spectrum and convergence relative to ibothros, respectively, 
for the turbulent accretion problem. This comparison is more challenging because the high 
energy emission originates in compact "hot spots." Evidently the two models agree at the 
10% level everywhere, with the largest differences at high and low frequency, where the 
subgrid reconstruction comes into play. Excluding these high and low frequency regions, the 
agreement between the codes is at the few percent level. 



8. Sample Calculation & Full Code Tests 

We now apply grmonty to the same HARM simulation data used in the grmonty- ibothros 
comparison above, this time with scattering enabled. Figure IT6l shows the resulting spectrum 
with the ibothros result shown for comparison. Evidently, scattering has little effect on the 
sub-mm spectrum since the scattering optical depth is small (~ 10~ 4 ) but the model now pre- 
dicts a significant X-ray flux. A more detailed analysis of Comptonized sp ectra from GRMHD 



simul ations in the context of Sgr A* will be given in a separate paper (IMoscibrodzka et al. 



20091 ) . 

Figure [17] shows the results of self- convergence tests for the full code. Convergence is 
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initially slow but quickly approaches a rate proportional to N~ 1 / 2 as spectral bins become 
sufficiently sampled. The spectra shown in Figure [161 which were taken from a single run 
with 10 9 superphotons, were used for references. 

9. Summary 

We have described and tested a code that solves the radiative transfer problem for opti- 
cally thin ionized plasmas in general spacetimes. The code treats the full angular dependence 
of emission and absorption, treats single Compton scattering exactly (double Compton and 
induced Compton scattering are neglected), and can be easily adapted to simulate emission 
from both analytic and numerical models. While we have specialized to synchrotron emission 
in this work, grmonty is constructed so that it is straightforward to include other relevant 
emission mechanisms such as bremsstrahlung with only minimal modification. 

As a demonstration of a practical use for grmonty, we have computed the first spectra, 
including synchrotron emission and Compton scattering, from GRMHD models of a turbu- 
lent accretion disk. Other potential applications of our code are to neutron star accretion, 
emission from relativistic blast waves, and any problem where relativistic bulk motion makes 
radiation transfer treatments that expand the flow in orders ofv/c problematic. 
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x [GM/c 2 ] 

Fig. 1. — Photon geodesies for isotropic emission from the rest frame of a fluid element in 
a marginally stable circular orbit around a Kerr black hole with a/M = 0.9375. Results 
shown from grmonty (points) and geokerr (lines). The point size varies linearly with the 
z-coordinate. 
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Fig. 2. — Average fractional error in the conserved quantities k t and as a function of step 
size parameter e. The solid line is e 2 , showing that grmonty's geodesic integrator converges 
at second order. 
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Fig. 3. — In the top panel, the spectrum of a synchrotron emitting sphere with low optical 
depth above ~ 10 8 Hz viewed nearly perpendicular to the magnetic field from grmonty 
(crosses) and from a semi- analytic procedure (solid line). The bottom panel shows the 
fractional difference between the two results. 
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Fig. 4. — Integrated fractional error in the grmonty spectrum for a synchrotron emitting 
sphere with low optical depth viewed nearly perpendicular to the magnetic field as a function 
of the number of superphotons produced. The results are similar for other magnetic field 
orientations. The dashed line is proportional to N^ 1 / 2 . 
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Fig. 5. — In the top panel, the spectrum of a synchrotron emitting sphere of high optical 
depth below ~ 10 11 Hz viewed nearly perpendicular to the magnetic field from grmonty 
(crosses) and for a semi-analytic procedure (solid line). The bottom panel shows the frac- 
tional difference between the two results. 
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Fig. 6. — Integrated fractional error in the grmonty spectrum of a synchrotron emitting 
sphere with high optical depth as a function of the number of superphotons produced. The 
dashed line is proportional to iV~ 1//2 . 
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Fig. 7. — Spectra (upper panel) from grmonty (points) and sphere (solid line) produced 
by Comptonization of soft photons in a homogeneous, spherical cloud of hot plasma. Com- 
putations are done for: plasma optical thickness r = 10~ 4 , plasma temperature O e =4 and 
the central source radiative temperature kT r /m e c = 10~ 8 . Lower panel shows the fractional 
difference between the two spectra. 
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Fig. 8. — Same as in Figure [7J but for r=0.1. 
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Fig. 9. — Same as in Figure [7J but for r=3.0. 
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Fig. 10. — Integrated fractional difference between grmonty and sphere for the spherical 
scattering test for optical depths of 10~ 4 (solid), 0.1 (short dash), and 3 (long dash). The 
dotted line shows the self-convergence results for the sphere code for an optical depth 
t = 10~ 4 . The dot-dash line is proportional to iV -1 / 2 . 
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Fig. 11. — Illustration of how interpolation can lead to a discrepancy between grmonty 
and ibothros when the spectrum is sensitive to grid-scale structure. Shown are the grid 
specified values for some fluid property (solid line), the interpolated values (dash line), and 
the average zone values based on interpolation (dotted line). 
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Fig. 12. — The top panel shows the spectrum of a radially infalling spherically symmetric 
source threaded with a radial magnetic field around a Schwarzschild black hole as computed 
by ibothros (solid line) and grmonty (crosses). The bottom panel shows the fractional 
difference between the two. 
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Fig. 13. — Integrated fractional error in the grmonty spectrum for the spherically symmetric 
Schwarzschild problem as a function of the number of superphotons produced. The dashed 
line is proportional to N~ l l 2 . 
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Fig. 14. — The top panel shows the spectrum of a snapshot from a HARM simulation of a 
turbulent accretion flow onto a Kerr black hole as computed by ibothros (solid line) and 
grmonty (crosses). The bottom panel shows the fractional difference between the two. 



-36- 




Fig. 15. — Integrated fractional error in the grmonty spectrum for the turbulent accretion 
problem as a function of the number of superphotons produced. The dashed line is propor- 
tional to N' 1 ' 2 . 
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Fig. 16. — Same as Figure [TU except Compton scattering is included. The histograms 
show the grmonty result for nearly edge-on and face-on inclinations and the solid line is the 
ibothros spectrum for a nearly edge-on inclination. 
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Fig. 17. — Self-convergence test results for the spectra shown in Figure fT6l Convergence for 
the edge-on (solid) and face-on (dot-dash) spectra are shown. The dashed line is proportional 



